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Abstract 

We present two improvements to the tight-binding approximation of time-dependent den¬ 
sity functional theory (TD-DFTB): Firstly, we add an exact Hartree-Fock exchange term, 
which is switched on at large distances, to the ground state Flamiltonian and similarly to 
the coupling matrix that enters the linear response equations for the calculation of excited 
electronic states. We show that the excitation energies of charge transfer states are improved 
relative to the standard approach without the long-range correction by testing the method on a 
set of molecules from the database in [/. Chem. Phys. 2008 , 128, 044118] that are known to 
exhibit problematic charge transfer states. The degree of spatial overlap between occupied and 
virtual orbitals indicates where TD-DFTB and lc-TD-DFTB can be expected to produce large 
errors. 

Secondly, we improve the calculation of oscillator strengths. The transition dipoles are 
obtained from Slater Koster files for the dipole matrix elements between valence orbitals. In 
particular, excitations localized on a single atom, which appear dark when using Mulliken 
transition charges, in this way acquire a more realistic oscillator strength. 

These extensions pave the way for using long-range corrected TD-DFTB (lc-TD-DFTB) 
to describe the electronic structure of large chromophoric polymers, where uncorrected TD- 
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DFTB fails to describe the high degree of conjugation and produces spurious low-lying charge 
transfer states. 


1 Introduction 

Tight-binding DFT (DFTB) 2 4 and its time-dependent formulation TD-DFTB 5 6 are semi-empirical 
methods based on (TD)-DFT 7 9 that inherit many of the advantages and shortcomings of the lat¬ 
ter. The failure of TD-DFT to describe charge transfer states 10 is particularly severe if one deals 
with extended molecules or oligomers with a large degree of conjugation. These charge transfer 
states, which appear at unphysically low energies, can be removed if a long-range exchange term 
is included, leading to the long-range corrected TD-DFT. 11 

In the tight-binding formulation this correction can also be included at almost no additional 
computational cost. The possibility to include a range separated functional into DFTB has been 
explored before by Niehaus and Della Sala. 12 They generate the pseudo atom basis starting from 
a full DFT calculation on a single atom with a range-separated functional. This changes the elec¬ 
tronic parametrization and makes the pseudoatoms depend on the long-range correction. Then 
they consistently add the first order correction arising from the long-range correction in the tight- 
binding Kohn-Sham equations and the linear response formalism. 

Here instead, a much simpler (maybe less rigorous) approach is proposed, where the existing 
electronic parameters (pseudo atoms and Hubbard parameters) are left untouched. Since the usual 
DFTB parametrizations are based on a local density approximation (LDA), the long range ex¬ 
change can be incorporated by simply adding the attenuated exact exchange energy E lr to the total 
electronic ground state energy. For the calculation of excited states, a long range correction term 
is added to the coupling matrix in the spirit of CAM-B3LYP. lj The exchange integrals are then 
approximated by products of transition charges as usual. The distance at which the exact exchange 
is gradually switched on is controlled by new a parameter R\ r (equivalent to 1 //i in the notation 
of ref. 11 ) that can be adjusted to fit excitation energies to CAM-B3LYP results or experimental 


2 


values. 


This article is structured as follows: In sections |2T| and |T2| the working formulae for the long- 
range correction and the transition dipoles are derived. Section |2.3| deals with different ways of 
quantifying charge transfer at the level of tight-binding DFT. Finally the method is applied to a test 
set of molecules and compared to CAM-B3LYP results in sections [3] and |4j 
Technical details are deferred to the appendices: 

• Appendix [A] elucidates ambiguities in the definition of the y-matrix that lead to differences 
between DFTB implementations. 

• Appendix [B] discusses different ways of splitting the Coulomb potential into a short and a 
long range part. 

• Slater-Koster rules for transition dipoles are derived in the appendices [C] [D] and |E| 

2 Method Description 

Many review articles have been written about DFTB,Q 4 16 for a pedagogical introduction see. 3 To 
make clear at which points we introduce modifications, we recapitulate here the basics of DFTB 
and TD-DFTB: 

To derive the DFTB Hamiltonian, one starts with the full DFT energy functional 16 

E\p] — ^fai^a | J ^ ext P^j l0a) + 2 J j +^xc [p] +£nuc (1) 

where (f> a are the Kohn-Sham orbitals of the system of non-interacting electrons with the occupa¬ 
tion numbers f a , V ext is the potential the nuclei exert on the electrons, as well as any external elec¬ 
tric field, E xc is the exchange-correlation functional and E nuc the energy due to Coulomb repulsion 
between nuclei. One proceeds by expanding the energy around a reference density, p — po + 8p, to 
second order in dp. The reference density po is a superposition of atomic densities of the individ¬ 
ually neutral atoms, while the redistribution of the electron density 8p is a result of the chemical 
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bonding. 


After rearranging the expansion of the energy functional, 


E[po + Sp] | + ^ext + Vh [Po] + V xc [Po]^ | 0a) 


1 

+ 2 



W + j^) Spip ' 


( 2 ) 


2 j ^n[po]po{r) + (^E xc [Po] - J V xc [po\po{r)) +E nuc , 


the energy is partitioned into terms depending on the reference density and the orbitals (SbsX the 
density fluctuations (£ C oui,xc) and the repulsive energy (S rep ) , which stands for everything else not 
covered by the first two terms: 


E[p 0 + Sp] «£/ a (0a | H[po] | 0 fl ) 

a 



+ E nuc + everything else 


band structure energy E/, s 

Coulomb + part of xc energy S cou l,xc 
repulsive energy S rep 

(3) 


To approximate S couPxc one assumes that the charge fluctuation 8p can be decomposed into 
spherically symmetric contributions centered on the atoms, 5p — Ag/F/(|r — Rj\), where Ac// 
are the excess Mulliken charges on atom /: 

W[P0 + Sp] = I/ /'(^I + _l_) w 

N at N at W 

= YY f ij 

/ . / , ■ C/ coul,xc 
/ J 

For partial charges sitting on different atoms, I ^ /, only the electrostatic interaction is taken into 
account (depending on the distance R/j of the atomic centers) and any exchange or correlation 
interaction is neglected, since the exchange correlation functional is assumed to be local. The 
interaction energy of charge with itself on the same atom is controlled by the Hubbard parameters 
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Uh, which can be obtained from experimental ionization energies and electron affinities, or ab- 
initio calculations. This leads to a partitioning of £ 77 ul xc ' nl ° picvcwisc contributions: 


pIJ 

^COUljXC 


jAq/Aqj f /' I^J 

\U H (A qi ) 2 I = J 


-AqiAqjyu(Ru) 


(5) 


The approximate DFTB energy is now a function of the partial Mulliken charges instead of the 
density: 


£dftb[{A<?/}] = | H[po] | (j) a ) + -^Yij(Rij)AqjAqj + {E[p] -£dftb)^ (6) 

a 1 i j s --- ' 

LkjVI^Ru) 

All the deviations from the true energy are absorbed into the repulsive potential, which ideally 
should only depend on the molecular geometry but not on the charge distribution. In a rough ap¬ 
proximation, these deviations can be decomposed into contributions from pairs of atoms, V^iR/j), 
and be adjusted to higher-level DFT methods. 

Fitting 17 and validating 18 the repulsive potentials is the most time-consuming part of parametriz¬ 
ing the DFTB method. Although properly adjusted repulsive potentials are crucial for molecular 
dynamics simulations, structure optimization or vibrational spectra, we will neglect them here, as 
they have no influence on the electronic absorption spectra. 

The Kohn-Sham orbitals are expanded into a minimal set of pseudo-atomic orbitals {| /i) }, 
which are compressed by a confining potential, since orbitals in a molecule are less diffuse than in 
the free atoms: 

« = E 41 #*) < 7 > 

A variation of the energy with respect to the Kohn-Sham orbitals, under the constraint that the 
orbitals are normalized, leads to the DFTB equivalent of the Kohn-Sham equations with the DFTB 
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Hamiltonian (in the basis of atomic orbitals p and v): 


N a , 


<v FTB = (M I H\po] I v) +- (p | v) £ (y IK + Yjk ) P e /, v g i 


( 8 ) 


■V*" 

A/0 


a:=i 




where p e / means that the atomic orbital p belongs to atom /. 

The Kohn-Sham equations are solved self-consistently: In each step the new partial Mulliken 
charges A cji and the new Hamiltonian are computed from the orbital coefficients of the previous 
iteration and the resulting Kohn-Sham equations are solved to give the next orbital coefficients. 
These steps are repeated until the charge distribution and the density matrix do not change any¬ 
more. 

Excited states are calculated in the framework of linear-response TD-DFT, which was adapted 
to tight-binding DFT by Niehaus. 5 For brevity, we only sketch the working equations of TD- 
DFT and how they are adapted to the tight-binding formulation. A converged DFT calculation 
provides the single-particle Kohn-Sham orbitals, and to a first approximation excitation energies 
are differences of virtual and occupied Kohn-Sham energies (0 OV(J = £ va — £oa- 

The coupling matrix 


K, 


ov<7,o v <y 


0o<r(n)0vcx(n) 


+ 


8 2 e x 


J?t-r 2 | 5p (7 (r 1 )5p a /(r 2 )_ 


ha' (h) ha' ir 2 )d i rid i r2 (9) 


represents the response of the Kohn-Sham potential to a perturbation of the electron density and is 
responsible for adding a correction to the Kohn-Sham orbital energy differences. 

In linear-response TD-DFT excited states are computed from the Hermitian eigenvalue prob¬ 
lem 8 

(A - B) 1/2 (A + B) (A - B ) 1/2 F l = Q.JF! (10) 
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where the matrices A and B contain the coupling matrix: 


A 


ova,of v 1 a' 


&o.o' < -V. v' 


&a,a' (£vc 


C-CKJ 


+ K n 


(ID 


B 


ova,o'v' o' 


= K, 


ova, v' o' a' 


( 12 ) 


The solution of eqn. (??) provides the excitation energies Ej = ft£t[ as eigenvalues and the coeffi¬ 
cients for single excitations from occupied to virtual orbitals (o —>• v) 


^ova - fQ 52 52 52 

v ' o'docc v'Gvirt a' 

as eigenvectors. 

In the absence of a non-local exchange term, (A — B) is effectively diagonal so that eqn. (??) 
can be simplified to yield Casida’s equation: 8 


(A-B) 1 / 2 


over .o r v r G r 


fL 


o'v'a' 


(13) 


^ ^ ^ &o.o'&v,v'$a.a' (Uct £oa) -(-2a/£ v(7 F^aK ova C) t \f£ v l ( y l £ 0 'a' 


o'aocc v'Evirt a' 


with the coefficients 


C = 

'-'ova 


£\ a £oa 

Cb 


Fo'v'a' — tfFova 
(14) 


(15) 


In the Tamm-Dancoff approximation 1920 excitation energies and coefficients are calculated 
from a different eigenvalue problem that results from eqn. (??) by setting B = 0: 


^ ^ £a ovct , oW c oVct , 


o Eoccv £virt a 


(16) 


Using the Tamm-Dancoff (TDA) approximation in conjunction with a long-range correction 
(which will be introduced later) can sometimes be advantageous for two reasons: 

• Firstly, it has been shown 21 22 that TDA excitation energies can actually be better than those 
obtained from the full solution of the LR-TD-DFT equation ??, in particular, when singlet- 
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triplet instabilities would lead to imaginary excitation energies. 


• Secondly, when (A — B) is not diagonal, the full solution requires the computation of the 
matrix square root (A — B) 1//2 which is computationally demanding unless one resorts to an 
iterative algorithm 23 specifically designed to deal with eqn. (??). 

The tight-binding approximation consists in replacing transition densities (j) 0 (r)(j) v (r) by transi¬ 
tion charges q°^ (defined later), so that the coupling matrix simplifies to 

Nat N a , 

Vv'^E (i7) 

A=\B=\ 

Similarly, the transition dipoles between Kohn-Sham orbitals are reduced to sums over transition 
charges on different atoms: 

(o | r | v) ss J^R A q^’ (18) 

A 

Note how in going from eqn. (??) to the eqn. (??) and in approximating the coupling matrix 
in eqn. (??) the exchange-correlation term has been neglected for charge distributions on different 
atoms, arguing that the xc-functional is local. This is where the long-range correction will be put 
to work. Eqn. (??) is usually a reasonable approximation, unless orbitals o and v are located on 
the same atom. 


2.1 Long-range correction for TD-DFTB 

The Coulomb potential is separated into a long-range and a short-range part, 24 25 where the posi¬ 
tion of the smooth transition between the two regimes is controlled by a parameter R\ r : 


1 

r 


| erf (ifc) 

r r 

' -V-'' V -V-- 

short range long range 


(19) 


The short range part of the exchange energy is treated with DFTB while for the long range part 
the exact Hartree-Fock exchange is used. Since in DFTB a local exchange correlation functional 


8 









is employed, the short range term is essentially neglected. 

The electron integrals of the screened Coulomb potential (for real-valued orbitals) 

/* C erf f —— 'j 

(ij\ab) b = I / 0/(ri)0 7 (n)— rid 3 r 2 (20) 

are approximated as in DFTB: 3 The transition densities between different orbitals p kl (j) — <j)k(7 )</>/(?) 
are decomposed into atom-centered contributions: 

p kl (j) = t P AW (2D 

A 

Next, the monopole approximation is made assuming that the transition density due to atom A is 
spherically symmetric around that center: 


p%(r)=q%F A (\r-R A \ ) (22) 

In fact, the exact form of the functions F A (r ) is not know. 3 They can be assumed to be Gaus- 
sianl 3 or Slater functions J 12 but in either case the width of the density profile should be inversely 
proportional to the chemical hardness (the Hubbard parameter U ) of the atom. 

With these approximations the long-range two-center integrals can be written in terms of tran¬ 
sition charges q } A and q c g\ 

r r erff^ 2 ) 

{ij\ab) Xx = YJL c ii^ / / ^a(|g ~Ra\) KRhJ F B (\r 2 -R B \)d 3 ri d 3 r 2 (23) 

A B J J G2 

The transition charges are defined as: 

Qa = \ £ £ fc c v + (24) 

z /ieA v 

where ,S' U v is the overlap matrix between the atomic orbitals /i and v and is the coefficient of 
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the atomic orbital /i in the molecular orbital i. 

Assuming a Gaussian- or a Slater-function form for F^{r) the integral for the unscreened 
Coulomb potential can be performed analytically 3 and the result is defined as the y-matrix: 

(25) 

./ ./ \n-r2\ 

For the long range part of the Coulomb potential the error-function can be taken out of the 
integral in eqn. ?? in a good approximation, giving the long-range y-matrix 

X\b( r ab ) « erf(^)y Afi (i? AB ) (26) 

/'ll 

from which the electron integrals with the long range part of the Coulomb potential can be calcu¬ 
lated as: 

(MA|<7v) lr ^££y^ A <7£ v (27) 

A B 

Incidentally, the authors of 12 have shown that the integral ?? can be transformed into a one¬ 
dimensional integral that can be solved numerically (assuming the /^(r)’s are Slater functions). 
As can be seen in Fig. [Tjthe approximate y^ g and the numerically exact solution differ in the limit 
Rab ~> 0 but coincide for larger distances. 

For atomic orbitals the transition densities are simply: 

E A) + E A))S^x (28) 

Here 5(/i e A) is equal to 1 if the atomic orbital /i is centered on the atom A and 0 otherwise. 
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Figure 1: Comparison between the exact solution of the integral in eqn. ?? and the approximation 
made in eqn. ?? for different (0 = The green curves are plotted for the long-range radius 
Rj, =3.0 bohr that is actually used in our parametrization. The approximation is quite accurate for 
interatomic distances that are relevant for charge transfer. For short atomic distances our approx¬ 
imation slightly overestimates and vanishes for R\k —* 0. In particular the diagonal elements 
of the Hamiltonian H^ v are not modified, since R,\k — 0 if the orbitals /i and v are located on the 
same atom. 


The long-range electron integrals in the basis of atomic orbitals can now be approximated as: 
(Mkv) lr « IZfog-qf = Wv»IE k* (%a) + «(UA))(J(<I6i) + S(vElt)) 

A B A B ^ 

= j$li?iSvc^^'Yab{8(L 1 e A)5 (ct g B ) + <5(/i g A)6(v g B) 

4 A B 

+<5(A g A)5(cr g B) + <5(A g A)5(v e B)} 

= 1 { r M<r + rjv + + rj v } 

(29) 


where the abbreviation T^ a = 'LaHb7ab^(I j - £A)S(<jeB) has been introduced. 

The additional contribution to the total energy E e i ec = F’dftb + due to the long-range ex- 
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change is: 


E *=-\ E P^PxAnMev)* 

4 ju,A,ct,v 


( 30 ) 


where /A, A, a and v enumerate atomic orbitals, and the density matrix elements / J ut7 are computed 
from the molecular orbital coefficients j as P^ a = 2^, eocc C u ,C tT( . 

Minimization of the total energy with respect to the density matrix leads to the Kohn-Sham 
Hamiltonian 


// 


KS 

MV 


r/DFTB i /r'lf 
n HV ' n IAV 


(31) 


with an additional term: 


// MV=-^^cr(iU^l^v)i r 

Z Aa 

=’ - f E^^m^o { r J«+ r !!v+rL+rj, 

Ac 

= -§|Ef r l:« (ewao )js- 

+ r li r vEE( S M f 'A ( ,)S«v 


+EE(v Km-Lijs 

cj A X V 77 

+ 52 *^/iA V.^Ag^crv^ I" 


(32) 


Av 


It is important to perform the sums over the indices in such an order that they can be implemented 
efficiently by nested matrix multiplications. 

The long-range contribution to the exchange energy in eqn. ?? can be computed as: 


E h — _ 

* 8 


§{e fE s M(E^vSv»)j^»rJ„+E (e s m*M (Ev s v„)r%,| (33) 


In the linear response formulation of TD-DFT the long-range correction leads to an additional 
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term in the coupling matrix, which shifts the excitation energies of the charge transfer states up. 
The corrections to the A- and B-matrices read (after separating the problem into separate singlet 
and triplet cases 2 ®): 


S B 


ov,oV 
ov,o'v' 


VM* &o) 4- 2 K ovo / v i + K qvo i v i 

2K ov yo> 4~ K ov,Vo! 


for singlets 


and 


^ov,oV VM* £o) T ^ov.o'v' 

T n _ L'\y 

D ov,o'v' — A ov,Vo' 

where the additional long-range coupling is given by: 


for triplets, 


(34) 


(35) 


<„V =~(°o>v , ) b ^-Z A LB‘f A °'rt B (.RAB)q'S' (36) 

=-K|vo')lr«-E4L««V/ B («4s)rf'’- (37) 


().()' are occupied and v, v' are unoccupied molecular orbitals. In addition to the transition 
charges between occupied and unoccupied orbitals, q™, which are needed for constructing the 
coupling matrix anyway, one has to calculate transition charges between occupied-occupied or¬ 
bitals, q °/, and between virtual-virtual orbitals q v g . 

In this approximation the quality of triplet excitation energies is expected to be much lower, as 
they are essentially equal to differences between Kohn-Sham orbital energies. Triplet states will be 
left aside in this work, since spin-unrestricted DFTB^ js necessary to describe them quantitatively 
and since they are dark in absorption spectra. 

The oscillator strengths of singlet states are obtained as 




4 

3 


I S> 

oEocc vZivirt 


2 


r | v) \Z&iC ] ov 


(38) 
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2.2 Slater-Koster rules for dipole matrix elements 

Transition dipole matrix elements determine the oscillator strengths in TD-DFT calculations of ex¬ 
cited states. In DFTB they are usually approximated by transition charges located on the individual 
atoms, 

{xif i \r\xir j ) = Y J R a q i L (39) 

a 

where R a is the position vector of atom a. 


V) 

= E4 

l 1 

1 A i(r-R^)) 

(40) 

w j ) 

= £4 

V 

| v(r-R v )) 

(41) 


are Kohn-Sham molecular orbitals in the basis of atom-centered numerical orbitals, and the transi¬ 
tion charges are defined as 


q l a = i L L ( c /i c v*V +c l v c J ^S v ^j (42) 

z /iea v v 7 

This approximation fails if a transition happens between molecular orbitals on the same atom, 
which are orthogonal by construction. In this case the overlap matrix simplifies to the identity 
matrix and the Mulliken approximation leads to vanishing transition charges. Therefore the excita¬ 
tion energies for localized n —> n* excitations are not shifted by the TD-DFT coupling matrix and 
no improvement of these energies relative to the Kohn-Sham energies is achieved. 27 Moreover, 
the oscillator strengths, which determine the shape of the spectrum, can sometimes be sensitive to 
the way transition dipole matrix elements are computed. For example, in reference^ 27 an on-site 
correction to the dipole matrix element is introduced, so that the 2 II states of nitric oxide, which 
would be dark using the Mulliken approximation, gain a small oscillator strength. 

However, the oscillator strengths can also be calculated without approximation from the tran¬ 
sition dipoles between atomic orbitals which are assembled using the Slater-Koster rules. The 
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technical details are explained in the appendix [C] (apparently this approach has already been im¬ 
plemented in other DFTB codes, but is not well documented in the literature). 

For transitions comprising orbitals on different atoms the oscillator strengths derived from the 
Mulliken charges and from the tabulated dipole matrix elements are very similar. However, when 
the transition is confined to a single atom, the Mulliken approximation can miss states which have 
weak oscillator strengths. As an example, consider the acrolein molecule, where the Si state is 
characterized by an excitation from the n to the n* orbital on the oxygen. Since the oxygen is 
much more electronegative than the carbon it is bound to, both the n and the n* orbitals consist 
largely of the atomic oxygen orbitals. The oscillator strength is very small, but it is not zero. Fig j2] 
compares the absorption spectra when the oscillator strengths are calculated from Mulliken charges 
or transition dipoles, respectively. 



Figure 2: Absorption spectrum of acrolein. The molecular structure was optimized at the 
BL YpEEE^ cc - p y X zP level, and a range-separation parameter of R\ r — 20.0 bohr was used in 
the lc-TD-DFTB calculation. With the Mulliken approximation (green curve) several states disap¬ 
pear, although they have a small oscillator strength. For example, the lowest excited states at « 3 
eV is of n —* n* character and has a tiny oscillator strength of 4 • 10 4 which cannot be seen in the 
green curve. 
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2.3 Diagnosing charge transfer states 

Visual inspection of the orbitals involved in an excitation is normally required to characterize an 
excited state as a charge transfer (CT) state. The authors of ref. 1 introduced a simple numeric 
test for detecting problematic excitations. They quantified the degree of spatial overlap between 
occupied and virtual orbitals of an excited state / using the quantity 

A= £ £ K) 2 /W)II'M?)|A (43) 

oeocc vevirt J 

where C ! ov is the coefficient for the single excitation from the Kohn-Sham Slater determinant that 
replaces the occupied orbital o by the virtual orbital v in the excited state /. 

They found that the error of DFT functionals without long-range correction such as PBE 31 and 
B3LYP,^ 32 correlates with the value of A, which can vary between 0.0 and 1.0. Values below 0.5 
indicate a charge transfer or Rydberg excitation, for which PBE and B3LYP will probably be in 
large error. 

Unfortunately, the integral over products of orbital moduli is difficult to calculate in DFTB 
without resorting to numerical integration. Therefore we replace the integral by 

O ov = j |0 o (r)| 2 |0 v (r)| 2 d 3 r (44) 


and define the new quantity 


a 2 = £ £ K) 2 

oeocc vevirt 


O n 


V O 00 0 v 


(45) 


that should behave similarly to A from eqn. ?? . Applying again the monopole approximation, 0 ov 
can be approximated by 


Oov 



F A (\r-R A \)F B (\r-R B \)d 3 r. 

-V”-^ 

&AB 


(46) 
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The overlap integral Q-ab of the spherical charge distributions centered on atoms A and B can be 
performed analytically assuming that F has a Gaussian profile: 


Ur) = 


(27TO-2) 


2^/2 


exp 


2cy i 


(47) 


where the width of the distribution is inversely proportional to the Hubbard parameter of atom A 
(see eqn. 29 of ref. 3 ): 


o A = 


1 1 

a /nU A 


(48) 


The integral is: 


Ub = 


(27r(o-2 + C jJ)) 


3/2 


exp 


1 1 


2 aj + ol 


Ra —Rb 


(49) 


For extended molecular systems the analysis of charge transfer by looking at the orbitals can 
become very cumbersome as the grids on which they are stored have to be very large. In the 
appendix [F] we provide an alternative analysis method in terms of the difference density between 
ground and excited stated, which is partitioned into atomic contributions. 


3 Computational Details 

The values for the confinement radius ro and the Hubbard parameter Uh that were used to parametrize 
the electronic part of DFTB are shown in Table [I] The Kohn-Sham equations with the PW92 33 
local exchange-correlation functional for the pseudo atoms were solved self-consistently using the 
shooting method. — In order to quickly assemble the DFTB Hamiltonian and the transition dipoles, 
Slater-Koster tables were generated for all pairwise combinations of valence orbitals for the atoms 
H,C,N and O. The parameter for the long range correction was set to R\, = 3.03 bohr (for compar¬ 
ison, CAM-B3LYPG2 uses R\ r — jr — 3.03 bohr, but only 65% exact exchange). So far we did not 
try to find an optimum value for R\ r . 

The structures of the test molecules were taken from the database 38 that has been published 
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Table 1: Electronic DFTB parameters. In general, the Hubbard parameter is set to the difference 
between experimental ionization energy and electron affinity, Uh — IE^EA, 36 where available, 
and the confinement radius is taken as 1.85x the covalent radius as reported in, 37 tq = 1.85r cov . 
Deviating from this rule, the Hubbard parameter of nitrogen was set to 0.530 Hartree (taken from 
the DFTB implementation of Hotbit ®), since there is no experimental value for the electron affinity 
available (the N anion is unstable). 


Atom 

ro / bohr 

Uh 1 Hartree 

H 

1.084 

0.472 

C 

2.657 

0.367 

N 

2.482 

0.530 

0 

2.307 

0.447 


by the authors of Ref. 1 The molecules N 2 , CO, H?CO and HC1 were excluded because the min¬ 
imal basis set of occupied valence orbitals used in DFTB is not suitable for describing Rydberg 
states. In the calculations labeled as lc-TD-DFTB the long-range correction as described above 
was included, Casida’s equation ?? was used for solving the linear response equations and the tran¬ 
sition dipoles were assembled from Slater-Foster files. In the calculations labeled as TD-DFTB 
the long-range correction was absent and the transition dipoles were calculated using transition 
charges. For comparison we also determined the excitation energies and oscillator strengths with 
full TD-DFT using the CAM-B3HYP 13 functional and the cc-pVTzJ 30 basis set as implemented in 
the Gaussian 09 3 ^ program package. 


4 Results 

Energies and oscillator strengths of specific excitations computed with TD-DFTB, lc-TD-DFTB 
and CAM-B3FYP are compared in Table[2j The HOMO-FUMO gaps are plotted in Fig. |4} 

Peptides. The first molecules in the test set are three model peptides of increasing length, 
a dipeptide, a /3-dipeptide and a tripeptide. The labels n(0,), 7r(A,) and nf, which are used to 
classify the excited states, refer to the lone electron pair on the the z-th oxygen, the binding K- 
orbital containing the z-th nitrogen and the anti-bonding ^-orbital on the z'-th carbonyl group. The 
ordering of the carbonyl groups and atoms is depicted in Fig. [3] 
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The energies of the lowest localized excitations are well reproduced both by TD-DFTB and lc- 
TD-DFTB. The energies of charge transfer excitations are underestimated, although lc-TD-DFTB 
puts them much closer to the CAM-B3LYP results and correctly predicts that the band of charge 
transfer states should be located above the local excitations. As opposed to this, TD-DFTB pro¬ 
duces charge transfer states that lie below the lowest local excitation. This problem gets worse as 
the peptide grows: In the tripeptide TD-DFTB erroneously predicts the lowest excitation to be a 
long-range charge transfer from n(0\) to 7r| with an energy of 4.59 eY, whereas the CAM-B3LYP 
energy of this excitation would be 8.68 eV. 



Figure 3: Peptide structures, a) dipeptide, b) /3-dipeptide and c) tripeptide. The carbonyl groups 
used to classify the excitations are encircled. 


Acenes. The next molecules are the smallest 4 polyacenes Naphthalene (n=l), Anthracene 
(n=2), Tetracene (n=3) and Pentacene (n=4). The large degree of conjugation leads to orbitals 
that are delocalized over the entire molecule. TD-DFTB underestimates the energies of all states 
consistently by < 0.5 eV, and this error remains stable with the size of the acenes. With lc-TD- 
DFTB the B 211 states deviate no more than 0.1 eV from the CAM-B3LYP reference values. 

N-phenylpyrrole. In this heterocyclic aromatic compound, TD-DFTB underestimates local 
excitations by « 1 eV while lc-TD-DFTB is correct to within 0.1 eV. Without long-range correction 
the lowest state with A\ symmetry has charge transfer character (A 2 = 0.02). The long-range 
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correction shifts this state to higher energies so that the lowest A\ state now belongs to a local 
excitation (A 2 = 0.49), as it should. 

The second Bi and A\ states, which involve an electron transfer from the pyrrole ring to the 
benzene ring, are predicted far too low in energy by TD-DFTB. 

DMABN. 4-(N,N-dimethylamino)benzonitrile (DMABN) possesses a low-lying charge trans¬ 
fer state that is formed when the nitrogen on one side of the phenyl ring donates charge to the -C=N 
group on the opposite side. Although the ordering of the states is correct even without long-range 
correction, lc-TD-DFTB comes much closer to the reference values than plain TD-DFTB. 

Polyacetylenes. As with the acenes, TD-DFTB energies are too low by 0.5 eV as compared to 
the lc-TD-DFTB and CAM-B3LYP values. 
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Figure 4: HOMO-LUMO gaps for the molecules from the test set. The names of the molecules 
encoded by the numbers on the x-axis can be found in the first column of table |2j As expected 
from a long-range corrected functional the HOMO-LUMO gaps are larger for lc-TD-DFTB than 
for TD-DFTB. 
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Table 2: TD-DFTB and lc-TD-DFTB excitation energies and reference values from TD-DFT cal¬ 
culations at the CAM-B3LYP/cc-pVTZ level. Energies are in eV and oscillator strengths are given 
in brackets. The type of excitation (L = local, CT = charge transfer, DL = delocalized) was as¬ 
signed by inspecting the dominant pair of occupied and virtual orbitals in the transitions. The A 2 
values are calculated for lc-TD-DFTB. 



Molecule 

State 

A 2 

Type 

TD-DLTB 

lc-TD-DLTB 

CAM-B3LYP 

1 

Dipeptide 

A” n(0\) -A K\ 

0.61 

L 

5.48 (0.000) 

5.33 (0.000) 

5.68 (0.001) 



A” n(0 2 ) -A 71% 

0.75 

L 

5.44 (0.000) 

5.40 (0.000) 

5.92 (0.001) 



A’ n(Ni) -A 7 T| 

0.08 

CT 

5.54 (0.002) 

6.23 (0.015) 

7.00 (0.010) 



A” n(0\) -A n\ 

0.20 

CT 

4.92 (0.000) 

5.98 (0.000) 

7.84 (0.000) 

2 

/3 -dipeptide 

A” n{0\) -A 7l\ 

0.77 

L 

5.47 (0.000) 

5.39 (0.000) 

5.67 (0.001) 



A” 11 ( 02 ) -A 

0.77 

L 

5.45 (0.000) 

5.39 (0.000) 

5.76 (0.000) 



A’ k(N\) -a 

0.56 

CT 

5.64 (0.000) 

7.37 (0.558) 

7.42 (0.328) 



A” n(0\) -A Ji% 

0.01 

CT 

5.05 (0.000) 

6.71 (0.000) 

8.38 (0.008) 

3 

Tripeptide 

A” n(0\) -A 

0.59 

L 

5.51 (0.000) 

5.34 (0.000) 

5.72 (0.001) 



A” n(0 2 ) -A ;r| 

0.58 

L 

5.50 (0.000) 

5.37 (0.000) 

5.93 (0.001) 



A” n(0 2 ) -A ^3 

0.72 

L 

5.44 (0.000) 

5.43 (0.000) 

6.00 (0.001) 



A’ 7T(Ai) -A 7T| 

0.07 

CT 

5.56 (0.002) 

6.25 (0.014) 

6.98 (0.014) 



A’ 7l(N 2 ) -A 7T 3 * 

0.16 

CT 

5.78 (0.002) 

6.52 (0.032) 

7.68 (0.102) 



A” n(Oi) -A 

0.19 

CT 

4.93 (0.000) 

5.99 (0.000) 

7.78 (0.000) 



A” n(0 2 ) -A ^3 

0.28 

CT 

5.16(0.000) 

6.33 (0.000) 

8.25 (0.000) 



A’ tt(Ai) -A 7T 3 * 

0.03 

CT 

5.20 (0.000) 

8.35 (0.000) 

8.51 (0.007) 



A” n(Oi) -A ^3 

0.01 

CT 

4.59 (0.000) 

9.04 (0.000) 

8.68 (0.000) 

4 

Acene (n=l) 

B 2 u 

0.83 

DL 

4.27 (0.007) 

4.53 (0.006) 

4.62 (0.000) 



B\ u 

0.77 

DL 

4.02 (0.044) 

4.84 (0.046) 

4.67 (0.071) 

5 

Acene (n=2) 

B\ u 

0.77 

DL 

3.00 (0.047) 

3.84 (0.067) 

3.53 (0.076) 



b 2u 

0.83 

DL 

3.66 (0.026) 

4.02 (0.021) 

4.04 (0.001) 

6 

Acene (n=3) 

B\u 

0.74 

DL 

2.32 (0.040) 

3.19(0.073) 

2.76 (0.071) 



b 2u 

0.82 

DL 

3.27 (0.056) 

3.69 (0.047) 

3.65 (0.003) 

7 

Acene (n=4) 

B\ u 

0.72 

DL 

1.85 (0.033) 

2.74 (0.076) 

2.22 (0.064) 



b 2u 

0.81 

DL 

3.01 (0.100) 

3.48 (0.087) 

3.39 (0.008) 

8 

N-phenylpyrrole 

b 2 

0.48 

L 

3.96 (0.005) 

4.98 (0.003) 

5.06 (0.013) 



Ai 

0.49 

L 

3.99 (0.000) 

5.07 (0.458) 

5.12(0.365) 



b 2 

0.29 

CT 

4.30 (0.009) 

5.24 (0.015) 

5.27 (0.015) 



Ai 

0.02 

CT 

4.51 (0.352) 

6.18(0.000) 

5.92 (0.179) 

9 

DM ABN 

B 

0.37 

L 

4.27 (0.023) 

4.47 (0.022) 

4.72 (0.024) 



A 

0.49 

CT 

4.52 (0.308) 

4.95 (0.453) 

4.91 (0.520) 

10 

Polyacetylene (n=2) 

B u 

0.68 

DL 

5.57 (0.500) 

6.21 (0.420) 

6.04 (0.706) 

11 

Polyacetylene (n=3) 

B u 

0.69 

DL 

4.54 (0.813) 

5.10(0.707) 

5.03 (1.110) 

12 

Polyacetylene (n=4) 

B u 

0.69 

DL 

3.88 (1.133) 

4.43 (0.995) 

4.39(1.533) 

13 

Polyacetylene (n=5) 

B u 

0.69 

DL 

3.41 (1.437) 

3.98 (1.269) 

3.94(1.961) 
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4.1 Correlation between errors and A 2 


When the deviations of the excitation energies relative to the CAM-B3LYP values are plotted 
against the degree of spatial overlap A 2 (see Fig. [5]), a clear correlation is visible. A 2 values close 
to 0.0 can be associated with charge transfer states, values around 0.5 with local excitations and 
values close to 1.0 with strongly delocalized excitations. For charge transfer states, the long range 
correction reduces the maximum error from -4.0 to -2.0 eV. For local excitations, TD-DFTB and 
lc-TD-DFTB have similar errors - TD-DFTB underestimates energies by at most 0.5 eV, while 
lc-TD-DFTB overestimates them by the same amount. For strongly delocalized excitations, as 
they occur in the acene series, lc-TD-DFTB shifts the error to the positive region with respect to 
TD-DFTB, and reduces somewhat the absolute values of the error. 


a) 


QJ 
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b) 


<u 

O 
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Figure 5: Deviation of excitation energies from CAM-B3LYP reference values for a) TD-DFTB 
and b) lc-TD-DFTB plotted against A 2 (a measure of spatial overlap defined in eqn. ??), local 
excitations (A), charge transfer excitations (•), delocalized excitations (■). The area into which 
errors from charge transfer states fall without long-range correction is highlighted in pink. 


Table [3] gives the mean errors averaged over the whole set of test molecules. In summary, 
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energies of localized, charge transfer as well as delocalized states are systematically improved 
by lc-TD-DFTB. In particular, the long range correction particularly markedly the energies of 
delocalized and charge transfer states while the energies of the localized states are only marginally 
better. 

Table 3: Mean errors (in eV) relative to the CAM-B3LYP excitation energies for the molecules 
in the test set. The long-range correction particularly improves energies of delocalized and charge 
transfer states. 

Method Total Local Charge Transfer Delocalized 
TD-DFTB LT3 051 2.22 046 

lc-TD-DFTB 0.46 0.34 0.83 0.17 

4.2 Are potential energy surfaces with lc-TD-DFTB accurate enough? 

The good agreement of lc-TD-DFTB absorption spectra with those obtained from full long-range 
corrected TD-DFT around equilibrium geometries unfortunately does not guarantee that global 
features of the excited potential energy surfaces are equally well reproduced. 

We have investigated the accuracy of the excited state potentials by exploring profiles of the 
excitation energies along the torsion angles in two of the molecules from the test set, DMABN 
and N-phenylpyrrole (see insets in Figs. [6] and [7] for the definition of the torsion angle t), and 
tried to identify the approximations in TD-DFTB that lead to the discrepancies with respect to full 
TD-DFT. 

Both molecules show dual fluorescence 4041 since the excited charge transfer state has an en¬ 
ergy minimum at a dihedral angle of 90°, while the ground state is more or less planar. As shown 
for DMABN, 42 density functionals without long-range correction such as PBE or B3LYP fail to 
predict an energy maximum of the charge transfer state at the twisted geometry. Only long-range 
corrected functionals are faithful to the shape of the potential energy surface. 

In order to separate the influence of the functional from the influence of the basis size, we first 
compare the torsion profiles for the PBE 31 with those of the LC-PBE functional 25 which is the 
long-range corrected version of PBE - both with the aug-cc-pVDZ 30 basis set (Figj6^) vs. b) and 
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Fig@i) vs. b)). In the next step we reduce the basis to a minimal one (STO-3G® 3 ) (Figj6j:) vs. 
d) and FigjTJc) vs. d)). Clearly, the absence of a long-range exchange term changes the profiles 
qualitatively, whereas the reduction of the basis size does not - except for a large energy increase 
of all excited states. Now, the stage is set for comparing with tight-binding DFT. 

In standard TD-DFTB, the charge transfer states collapse to a minimum at z — 90° for both 
molecules as expected of a functional lacking long-range exchange (Figj6jf) and Figj7p). Density 
functionals with local exchange cannot describe charge transfer over long distances because the 
transition density 0 o (r)0 v (r) vanishes if the occupied and virtual orbital are spatially separated. 44 

For lc-TD-DFTB, one would expect an improved description of the charge transfer states. 
While the energy maximum at the twisted geometry is at least qualitatively reproduced for N- 
phenylpyrrole (Figj7|e),g)), for DMABN lc-TD-DFTB performs even worse than TD-DFTB (Figj6]e)): 
The lc-TD-DFTB torsion profile shows an even more pronounced dip at 90°, but for a different rea¬ 
son. This becomes evident from Tables [4] and [5] which show the blocks of the long-range part of 
the coupling matrix Kj r that mixes different transitions between frontier orbitals. The lowest two 
excited states in DMABN are dominated by excitations from the HOMO-1 or HOMO into the 
LUMO or LUMO+1. A small off-diagonal matrix element of —0.01 is present at z = 0°, whereas 
at z — 90° the block becomes almost diagonal (the largest off-diagonal element is below smaller 
than 0.0001). Therefore at the twisted geometry, the coupling matrix cannot correct the excitation 
energies. The vanishing of the long-range coupling is probably due to the Mulliken approximation 
that reduces the transition density to spherically symmetric charge distributions around atoms. 

Table 4: long-range coupling matrix Kj r for the frontier orbitals at z = 0° for DMABN. 


K h . 

H-l -+ L 

H —)■ L 

H —> L+l 

H-l —>■ L 

0.1798 

0.0000 

-0.0106 

H —)■ L 

0.0000 

0.1543 

- 0.0000 

H —)■ L+l 

-0.0106 

- 0.0000 

0.1721 


Comparing the torsion profile of N-phenylpyrrole for the aug-cc-pVDZ and the minimal STO- 
3G basis sets allows one to identify a second source of error. Diffuse orbitals probably play a large 
role in charge transfer states, where orbitals participating in a transition overlap only slightly with 
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Table 5: long-range coupling matrix K[ r at T — 90° for DMABN. 


Kjr 

H-l -> F 

H —)• F 

H —>F+l 

H-l —> F 

0.1520 

-0.0001 

-0.0001 

H —)• F 

-0.0001 

0.1590 

- 0.0000 

H —>■ F+l 

-0.0001 

- 0.0000 

0.1690 


each other’s fuzzy tails. The torsion profile at the LC-PBE/STO-3G level (Fig. [7]c)) is comparable 
to lc-TD-DFTB (Fig. [7]e) and d)), which also uses a minimal set of valence orbitals. In both cases 
the maxima of the lowest 3 excited states are less pronounced than in the FC-PBE/aug-cc-pVDZ 
(Fig. |7]a)) case. 
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Figure 6: Excitation energy profile of DMABN along torsion angle. a),b),c) and d) with full 
TD-DFT using the PBE functional; e) and f) with tight-binding DFT. 
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Figure 7: Excitation energy profile of N-phenylpyrrole along torsion angle. a),b),c) and d) with 
full TD-DFT using the PBE function; e),f) and g) with tight-binding DFT. In e) the lc-TD-DFTB 
spectrum has been obtained from the full solution of eqn. ?? while in g) the Tamm-Dancoff 
approximation was employed. 
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5 Conclusion and Outlook 


Tight-binding DFT is a very efficient method for calculating the electronic structure of large or¬ 
ganic molecules around their equilibrium structures. This efficiency rests on the precalculation of 
matrix elements using the Slater-Koster rules and on approximations for the two-center integrals 
using Mulliken charges. Any improvement of the method should not deviate from these principles 
as it would forfeit speed of execution. We have presented two improvements that add very little 
computational overhead, the long-range correction in the spirit of CAM-B3LYP and the calcu¬ 
lation of oscillator strengths from Slater-Koster files for transition dipoles. In full TD-DFT, the 
computational cost increases considerably if range-separated functionals are used, not so in the 
tight-binding approximation. 

Two measures of charge transfer were tested which can be extracted easily from TD-DFTB 
calculations, the quantity A? and the particle-hole separation. These can be used to detect cases of 
problematic behaviour involving long charge-transfer character of excited states. Application of 
the method to a test set of molecules with problematic charge transfer states led to an average error 
of 0.3 eV for local excitations, 0.8 eV for charge transfer excitations and 0.2 eV for delocalized 
excitations. Especially energies of charge transfer states are improved by over 1.0 eV relative 
to tight-binding DFT without range-separation. A more general test of the method for a set of 
common organic molecules 45 is compiled in the supporting information. 

Of course, being a semi-empirical method, the errors of lc-TD-DFTB also depend on the right 
choice of transferable parameters and more work needs to focus on optimizing these parameters, 
in particular the range-separation radius R\ r and the choice of the switching function (see appendix 
|B} For simplicity, we derived the electronic parameters (the confinement radius ro and the Hubbard 
parameter Uh) from experimental data. Previous work on the parametrization of DFTB 17 46 47 has 
proven that, by fitting the parameters to full DFT calculations on training sets of representative 
molecules, the accuracy of ab-initio electronic structure calculations can be reached. 

In the future, we intend to use lc-TD-DFTB to study polymers containing e. g. organic chro- 


28 




mophores. In this context, a comparison with other semi-empirical methods is pertinent. While 
some approximations are similar to semi-empirical wavefunction based methods, the derivation 
from time-dependent density functional theory is advantageous from the point of view of compu¬ 
tational efficiency since only single excitations are included. This approximation is often sufficient 
for describing low-lying bright excited states. The problem that the active space grows quickly 
with the system size becomes relevant only for much larger systems, when the coupling matrix 
cannot be diagonalized with all single excitations included. Therefore the method is suited for 
predicting absorption spectra of large systems with a few thousand second row atoms. 
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A Different approximations for the y AB matrix 

Reproducing DFTB results requires some care not only due to the different parametrizations but 
also because the various implementations employ slightly different approximations. One such 
ambiguity concerns the definition of the y-matrix: 

Tab = [ f F A (\r l -R A \)—^—F B (\r 2 ~R B \)d 3 r l d 3 r 2 (50) 

J J \n-r 2 \ 

where R = \R A —R B \. 

The y-matrix describes the change of the Coulomb energy due to charge redistribution between 
the atoms A and B. Charge fluctuations are assumed to be spherically symmetric around each atom 
but the exact functional form is unknown. Slater and Gaussian functions are obvious candidates, 
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since the Coulomb integrals are well-known for these functions. In some DFTB-implementations, 
the y-matrix is based on Slater functions (presumably DFTB+ and the older code by Seifert) while 
in others it is based on Gaussian functions (HotbifS). 


A.l Slater functions 

If the charge fluctuation around an atomic center A is modelled by a Slater function, 

F A (\r-R A \) = ^exp (-T A |r-i? A |J (51) 

the Coulomb integral between two such charge distributions separated by a distance R = R A —Rb\ 
reads: 


Jab 


x a (x 2 {2 + t a R)-x 2 (6 + x a R)) _ XaR 

~ R[ 2(rj-rj) 3 

_ x a a {xI(2 + t b R)-x 2 b {6 + x b R)) _ XbR 

2 (^!) 3 £ - 


(52) 


In the limit, that the charge distributions are centered at the same position, the y-matrix becomes: 


lim Yak = 
R-+ o 


x A x B (xl + 3x A x B + xl) 


2(Ti + Tj)^ 


The width parameter t a is fixed by the requirement y AA (R — 0) = U A : 


(53) 


16 

Xa = y Ua 


(54) 


If both atoms are of the same type, one finds: 


lim Yab(R) 


1 

R 


4S + 33{xR) + 9{xR) 2 + {zR) 3 _ zR 
48 ^ 


(55) 
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A.2 Gaussian functions 


In the case of Gaussian functions, the charge distribution is modelled by 

^ (l ^ l) = (^ exp p^) <56) 

Again, the Coulomb integral between two Gaussian charge fluctuations can be performed analyti¬ 
cally, yielding: 


with 


Oa is determined by 


Jab 


Cab 


erf (C ab R) 



R 


1 


(Oa + °b) 


lim Jaa 
r^o 


lim 

/?->() 



R 


1 

Vkoa 


U a 


(57) 


(58) 


(59) 


which leads to 


1 

aA ~7m, 


(60) 


for the width parameter. 


B Long-range switching functions and the R AB —> 0 limit 

Our long-range correction for DFTB has been criticized for neglecting the short range correction to 
the exchange energy and using an approximate Yab that has the supposedly wrong R A b —» 0 limit. 
Both issues are related and we would like to defend our method by the following argument: 

The error function is not the only and probably not the best way to separate the Coulomb 
potential into a short and a long range part. For a general switching function /(r) the separation 
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becomes: 


1 

r 


1 -m t m 



short range long range 


( 61 ) 


Iikura’s LC-functionals 25 and Niehaus’ range-separated DFTB employ the error function 


/erf (r) = erf (cur) 


with CO — —. 

/?lr 


(62) 


But there are other choices. Toulouse^ has shown that adding a Gaussian to the switching 
function (termed erfgau ) 


2 ( 1 o 

/erfgau (r) = erf (cor) - -^=(cor)ex p (--((Or) 


(63) 


has the advantage that the correction to the exchange functional that needs to be added to compen¬ 
sate the presence of the long-range HF exchange at short distance is reduced. An expansion of the 
exchange functional around (o = 0, 


£"erfg a „ = £, + 0(ffl 5 ). (64) 

demonstrates that the correction is negligible for co < 1. 

The limit of the long-range y-matrix for Rab —» 0 depends also on the switching function. We 
approximate this matrix as: 

y/ s = J j F A {\n-RA\) f -^F B {\h~RB\)d 3 nd\ 2 

« f(R AB )J y > F A (|r 1 -^|)-|-F fi (|r 2 - J R fi |)J 3 r 1 J 3 r 2 (65) 

= f{RAB)jAB 

One can argue that the approximation is not warranted, since for erf the numerical solution of 
the integral leads to a maximum for Rab —> 0 (see Fig. 1 in reference ^ while our approximation 
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gives 0. However, if erfgau is used as a switching function the integral approaches a minimum in 
the limit Rab 0. The shape of ^'(Rab) at small distance depends strongly on the choice of the 
switching function, while the large distance limit is not affected. 

We proceed by computing /'" numerically for the erf and erfgau switching functions and com¬ 
pare the numerically exact results with the approximation in eqn. ??. It should become clear that 
the short range behaviour is immaterial and that it is actually beneficial that the long-range correc¬ 
tion vanishes for Rab —> 0, as then the short-range modifications to the exchange functional (or the 
error incurred by neglecting them) are smaller. 

Assuming a spherical distribution 

2 

F a = j^exp (66) 


the integral can be transformed into a one-dimensional integral (see Niehaus 12 ) 


-7-4^4 i 
„4r _ Z A Z B 1 
1AB — —F, 7 

ttRab 1 


sin (q R ab ) 


lo (q 2 + T 2 ) (<? 2 + T 2 ) 


J(q)dq 


(67) 


where f(q) is the Fourier transform of the switching function: 



f(r)e lqr dr 


assuming f(-r ) = -f(r) 


( 68 ) 


The Fourier transforms of the two switching functions are 

/«*(?) = P (-^) (69) 

and 

2i / q 2 \ 3 3 / 2 i / q 2 \ 

/„, s »(9) = -“P (-4^J - (-3^ ) P0) 

In Fig. [8]the switching functions and /'"-integrals are compared. The short distance behaviour 
is completely different for erf and erfgau. 
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erfgau 





Figure 8: In a) and b) the separation of the Coulomb potential into a short and long range part is 
depicted for the switching functions erf and erfgau, respectively, c) and d) show the resulting dis¬ 
tance dependence of j r for the two switching functions. The dashed lines show the approximation 
according to eqn. ?? and the solid lines show the numerically exact integrals ??. The thin lines in 
d) show the erf-integrals for comparison. 


C Precalculated dipole matrix elements 

The dipole matrix element between two atomic orbitals /i and v can be rewritten, so that it depends 
on the relative position R^ v — R v —R^ of the two centers and the overlap of the two orbitals: 

(H(r-Rn) \r | v(r-R v )) = Jd 3 r(j)*(r-R^)r^) v (r-R v ) 

= d 3 r'(j)*(r) ^+ R^ +R^-R V ) (71) 

= Jdr 3 ^(r)ty v (r -R^ v )) 

The second summand in the last line of eqn. ?? can be calculated from Slater-Koster rules for 
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the overlap matrix elements between orbitals /i and v 


R 





( 72 ) 


whereas the first summand needs special treatment as it contains the dipole operator. The Cartesian 
components of the dipole operator can be written in terms of p x , p y and p z orbitals located at the 
origin: 



(73) 

(74) 

(75) 


Therefore the dipole matrix element can be understood as the overlap of three orbitals at two 
centers: 

• the orbital /i at the origin, 

• a vector of p-orbitals representing the direction of the position operator, also centered at the 
origin 

• and the orbital v at the center v = R v — R^. 

The atomic orbitals consist of a radial part R„i(r) and an angular part T/ jW (0,0): 


u(^) — 0) 


(76) 


The angular parts are real spherical harmonics: 
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Orbital 

l 

m 

Yim(e^) 

s 

0 

0 


1 

20r 

Py 

1 

-1 

1 

2 i 

/ |sin(0)sin(0) 

Pz 

1 

0 


2 V / J COS ( 0 ) 

Px 

1 

1 

1 

2) 

/|sin(0)cos(0) 


The radial part R n i(r) is specific to each atom type and is obtained by numerically solving 
the radial Schrodinger equation for the atomic Kohn-Sham Hamiltonian with a local exchange 
correlation functional. For the valence s- and p x ,p y and p z orbitals ( l,m ) would take the values 
(0,0), (1,1), (1,-1) and (1,0) respectively. The valence shell of carbon, for instance, requires 
two radial functions, R^ =2 /=0 (r) for the 2s orbital, and R^ =1 / =1 (r) for the three 2p orbitals. 

Decomposing the orbitals /i and v into their radial and angular parts and expressing the dipole 
operator by eqns 73 174| and [75] the integral for the first summand in eqn. ?? becomes (in spherical 


coordinates): 


• poo pH p2K 

d 3 r f ru(r)7'Mr-Rny)= / r\dr\ / sin(0i)<i0i / d(j)\ 


i o 


1 0 


1 0 


d^r\ 


/ 4 jc 

x ,(n)F, <h)d T n 

^—-V-' V J 

(o) 


^ Y lA (d h (j)i) ^ 
?l,-l(01,0t) 


Rn v I v ( In -Rnv\)Yl v m v (^2:<h) 

^ --v---^ 

<t>v(r 2 ) 


n 


(77) 


where r 2 — \ r\ — Ru v \ ■ 02 and (f >2 depend on the integration variables n, 0 1 and (j)\ as illustrated in 
Figj9ji). 

We wish to find a way to precalculate these integrals and tabulate them, so that at runtime no 
integrals have to be solved. At first it seems, as if one had to solve the integral for all possible rel¬ 
ative arrangement in 3D space (expressed by R^ v ) of the two orbitals. Slater-Koster rules 48 allow 
to break the integral down to a set of few elementary integrals, that only depend on the relative 
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distance, and from which the integrals for any relative orientation can be assembled quickly. 

To derive the Slater-Koster rules for dipole matrix elements, we begin by rotating the coordinate 
system such that R 2lv points along the z axis. Since spherical harmonics form a representation 
of the rotation group SO( 3) for each angular momentum /, the action of this rotation is to mix 
spherical harmonics with different m but the same /. For (complex) spherical harmonics the mixing 
is described by the Wigner D-matrices, analogously real spherical harmonics, as they are used 
for orbitals, will be transformed by combinations of those ZTmatrices, which will be called D- 
matrices: 

i 

■<!>)]= £ (78) 

m!=—\ 

The D matrices are expressed as functions of three Euler angles A, B and F. In order to align the 
vector R Ll v (whose spherical coordinates are R. 0 and <f>) with the z-axis, the angles have to be set 
to A — j, B — 0 and T = <f>. In the integral ?? all three spherical harmonics have to be rotated 
leading to: 


d i r^{r')r<^ v {f-R^ v )= £ £ £ 

m\=—lfx 1112=—l WI3 =—Z v 


nifi ,tn 1 


( d\ ^ 

1 ,m 2 

\ [) ln 2 ) 


D lv 

m v ,m2 


T i{ ln ,ly ,my i m \ i m 2 > m 3) 


X Jo ridn Jo Sin (^)^ 1 ^^(n)^ V , / v(^ 2 ) 

/ 4 tz r~ K 


(79) 




The rotated coordinate systems and mixing of spherical harmonics is depicted in Figj9jt) and c). 

Since after the rotation the z-axes for both orbital centered coordinate systems coincide, one has 
( pi — (j>2 and the 0 integral can be done analytically. 49 For s and p orbitals 8 different expressions 
1 ,^ 2 ) result, which are listed in the appendix |D| in table [6] The remaining two-dimensional 
integrals over r\ and Q\ are best performed in cylindrical coordinates p,z. The variable transfor- 
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mations are given by 


r\ 

= \/p 2 + (z-h) 2 

(80) 

r 2 

= \Jp 2 + (z + h) 2 

(81) 

0i 

— arctan2 (p,z — h) 

(82) 

02 

— arctan2 (p,z + h) 

(83) 


where h = Rj ^- is half the distance between the two atomic centers (see Figjojl) ). Since the orbitals 
are highly peaked at the atomic position and decay exponentially towards larger distances, a grid 
with sampling points clustered around the two atomic centers (see Fig. B> ) is most suited for 
accurate quadrature. 

In total, one has to calculate the 8 two-dimensional integrals Dj(R^ v ) listed in the third column 
of table [6] as a function of the orbital separation and save them to a file. The coefficients Tj(x,y,z ) 
depend on the directional cosines of the vector R^ v , i.e. x = cos(a),y = cos(/3) and z — cos(y) as 
labelled in Fig j9}i), and account for the relative orientation of the orbitals. 


d 3 r / <j>*(r f )r f <l>v(r f —R[iv) ^Y^Ti(x,y,z)Di{R 




(84) 


The Slater-Koster rules for computing 7,0, are given in table [7] of the appendix [E| The 3 
integrals needed for dipole matrix elements between the valence orbitals of carbon and hydrogen 


are shown in Fig. 10 as a function of the interatomic distance. 
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Figure 9: a) Local coordinate systems 1 and 2 around orbitals /i and v. The spherical coordinates 
of a position r' can be specified with respect to either axis. The direction of the vector R^ v joining 
the two orbital centers is defined by the directional cosines x = cos (a), y = cos(/3) and z — cos(y). 
b) After rotating the coordinate systems R^ v coincides with the z-axes. c) The rotated spherical 
harmonics are linear combinations of spherical harmonics aligned with the axes, d) Cylindrical 
coordinates, e) Grid for integration. Two polar grids centered at the atomic positions are merged 
for an efficient distribution of sampling points around both atoms (implementation of Hotbit ®). 
The plane is divided into rectangles or triangles, whose size increases away from each center. 
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Figure 10: Precalculated dipole integrals for the atom pair h-c that are tabulated in Slater-Koster 
files, D\ (blue), Z) 2 (red) and D 2 (violet). 


D Unique radial integrals D t 


Table 6: List of the angular and radial integrals. /?„ / (r) denotes the radial function with angular 
momentum for the orbital p of one atom type. ri,r 2 and 0\ and 0 2 depend on the cylindrical 
integration variables p and z as explained in eqns. 


i 

0/(01,02) 

Integrals A(Aiv) 

1 

^sin(0|)sin(0 2 ) 

f j dzpdpR,^ ,/ M=0 (n )rii?,, Vl / v =i (r 2 )0 1 (0i, 0 2 ) 

2 

-fd c °s(0l) 

//rfzprfpi?,i M5 / M=o (ri)ni? Mv; / v=o (r 2 )0i(0i,0 2 ) 

3 

4 ^cos( 0 i)cos( 0 2 ) 

f f dzp dpR n ^ = o (n ) r\R nvi i v= \ (r 2 ) 0i ( 0i, 0 2 ) 

4 

87^ sin2 ( 0 i) 

//rfzprfpi?„ M; /^i(ri)ni? Mv;/v = o (r 2 )0i(0i,0 2 ) 

5 

I\ZI sin2 ( 0 i) cos (0 2 ) 

S SdzpdpR n ^i(n)r\R nv t v= i(r 2 )$\(e\, 62 ) 

6 

| \J\ sin(0i)cos(0i) sin(0 2 ) 

f fdzpdpR „ M ,/ M= i (ri)riR„ Vjly=l (r 2 )</> 1 (0i, 0 2 ) 

7 

jTjcos'ce,) 

J/Jzpc/pi?„ M5 / (ll= i(r 1 )rii?„ V) / v =o(r 2 )0i(0i,0 2 ) 

8 

It|cos 2 (0i)cos(0 2 ) 

//rfzprfpi?^ !/M= i(ri)ni? Mv)/v =i(r 2 )0i(0i,0 2 ) 


80-81 
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E Slater-Koster rules for assembling dipole matrix elements 


Table [7] summarizes the rules for calculating dipole matrix elements between two atomic-centered 
orbitals from Slater-Koster tables for the integrals D l (R\ 2 ). The atomic orbital <j)\ is centered on 
an atom at position R\ and the second atomic orbital </b is centered on another atom at position AS- 
x.y and z are the directional cosines for the vector R \2 pointing from the first to the second center. 
The 8 integrals A'=t,...,8( r ) arc precalculated for all distances r — |/?] 2 1 and tabulated. So far only 
s- and p-orbital are considered. Since the rules were obtained from a computer algebra program 
written for the software package Mathematical 0 rules for ri-orbitals could also be obtained easily. 
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Table 7: Slater-Koster rules for dipole matrix elements (for s and p-orbitals). For clarity the dis¬ 
tance between orbital centers is named r instead of v . 


orb. 

coord. 

orb. 02 

Rule for (0i(r-^i) | (r-Rn^ \ <h.(r-R 2 )) 

s 

y 

s 

yDi{r) 

s 

y 

Py 

(. x 2 + Z 2 )D l (r)+y 2 D 3 (r ) 

s 

y 

Pz 

yz{D 3 (r) -Di(r)) 

s 

y 

Px 

xy(D 3 (r) — D\ (r)) 

s 

z 

s 

zD 2 (r ) 

s 

z 

Py 

yz{D 3 (r)-D x (r)) 

s 

z 

Pz 

(x 2 +y 2 )Di(r) +z 2 D 3 (r) 

s 

z 

Px 

xz(D 3 (r) —D\(r)) 

s 

X 

s 

xD 2 (r) 

s 

X 

Py 

xy(D 3 (r) —D\(r)) 

s 

X 

Pz 

xz(D 3 (r) —D\(r)) 

s 

X 

Px 

(y 2 +z 2 )D l (r)+x 2 D 3 {r) 

Py 

y 

s 

{x z +y z )D 4 (r)+y 2 D 1 (r) 

Py 

y 

Py 

y{x 2 +z?)(D 5 (r)+2D 6 (r))+y 3 D s (r) 

Py 

y 

Pz 

z((x 2 +z 2 )D 5 (r)+y 2 (D 3 {r)-2D 6 {r))) 

Py 

y 

Px 

x({x 2 +z 2 )D 5 {r)+y 2 (D s {r)-2D 6 (r))) 

Py 

z 

s 
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Py 

z 

Py 
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Py 

z 

Pz 
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Py 

z 

Px 
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Py 

X 
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Py 

X 

Py 
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Py 

X 
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Py 

X 
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Pz 
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Pz 

z 
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Pz 

z 
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Pz 

z 

Pz 
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Pz 
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F Analysing charge transfer with density differences 


Different “metrics” have been devised to identify charge-transfer states automatically and warn 
about possible failures of TD-DFT. Since the A-metric 1 cannot detect all problematic charge trans¬ 
fer excitations, Guido et.al. introduced the Ar 51 - and r 52 -metrics. In particular the Ar-metric has 
an intuitive interpretation as the electron-hole distance. 

We can define a similar quantity at the tight-binding level. To this end, we start with the density 
difference between the excited state and the ground state, Ap/ = P/ — p 0 . In the linear response 
regime the Kohn-Sham “wavefunction” of the excited state I is a linear combination of single 
excitations from the Kohn-Sham ground state Slater determinant: 

I Vi) = £ I Ci 0 ala 0 | T'o) (85) 

vGvirtoGocc 

The operator for the electron density in second quantization reads (in the basis of Kohn-Sham 
orbitals): 

P (?) = £ £0a (?) 0/3 (?) (86) 

a p 

Here, o and o' denote occupied, v and v' virtual and a and /3 general Kohn-Sham orbitals. 
Combining eqns. ?? and ??, the density of state I becomes: 

Pl(r) = {'i'j|p|f/) = ££E «(?% (?)c£c; v ('*' o | alMUpaU* I fo) (87) 

o,o' v,V a,p 

By using anti-commutation relations for Fermions, and the fact that the ground state only contains 
occupied orbitals (so a v \ 'To) = a ' 0 \ 'To) = 0), the expression for the density can be reduced to: 53 

Pl(f) = EE c »4„«?)'C..(7)-EEc!:c; < ,«(r)fc(r) + £|^(?)| 2 (88) 

O v,v' V 0,0 1 o 

= Pp (?) + P/j (?) + Po (?) (89) 
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In the exciton picture, the /-th excited state can be described by a bound particle-hole pair. The first 


term in eqn. 88 can be identified with the particle, the second term belongs to the hole, while the 
last term is just the ground state density. The difference density between the ground and excited 
state is the sum of the particle and hole densities: 


Ap I = p,-p 0 = p e + p h (90) 

Since tight-binding DFT deals with partial charges instead of a continuous density distribution, 
we have to coarse grain the density to an atomic resolution. Again, the transition charges are 
approximated as a sum of spherically symmetric charge distributions centered on the individual 
atoms: 


0v(r)M?) = F a{V-Ra\) 

A 

(91) 

€(r)Mr) = ^FAV-Ra]) 

(92) 


A 


The particle and hole densities are partitioned into atomic contributions, 


Pe = iEc::cU'A(i-^i)=i^(i-^i) (93) 

A o.y.v' A 

Pl, = Y,Y,Cl:ci 0 ,q%’’F A (\r-R A \) = Y 1 <l h A FA(\?-RAl) (94) 

A v,o,o ' A 


q A and q h A are the particle-hole charges. 

The average positions of the particle and the hole result from the weighted average of all 
charges. 


h 


LaVa^a 

LaVa^a 

Uq h A 


(95) 

(96) 
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The particle-hole separation <i e _/, = | r e — r h \ indicates the spatial extent of the exciton. FigjTT] 
shows that in the /3-dipeptide, absorption of a photon can lead to a separation of positive and 
negative charge over a distance of 10.0 bohr. 



Ap 


Figure 11 : Difference density Ap for the 3rd excited state of the /3 -dipeptide. An exciton is formed 
as an electron jumps from the first carbonyl group to the second leaving a positive hole behind. 
The radii of the red (blue) circles are proportional to the hole (particle) charges on each atom. 


G Supporting Information 

A more extensive assessment of the quality of lc-TD-DFTB excitation energies is given in the Sup¬ 
porting Information. This material is available free of charge via the Internet at http://pubs.acs.org 
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